% kpmf using gradient descent
function [U, V] = kpmf_grad_descent(R, mask, D, K_u_inv, K_v_inv, sigma_r, eta, maxIters, minIters, thresPercent)

N = size(R, 1);     % #rows
M = size(R, 2);     % #columns

denomCompU = K_u_inv + K_u_inv';
denomCompV = K_v_inv + K_v_inv';

%**************************************************************************

% Initialize U and V

U = 0.1 * rand(N, D);
V = 0.1 * rand(M, D);

% FIXME: seems no need to use kernels in initialization
% U = zeros(N, D);
% V = zeros(M, D);

% for d = 1 : D
%     U(:, d) = mvnrnd(zeros(1,N), K_u); 
%     V(:, d) = mvnrnd(zeros(1,M), K_v);
% end

%**************************************************************************

% Start learning
numIters = 0;
stopUpdateAtIter1 = 0;
stopUpdateAtIter2 = 0;
stopTimes = 0;
while(1)
    
    % Gradient descent on U
    U_new = zeros(N, D);
    
    UV=U*V';  
    pDeriv=eye(N);
    for d=1:D
        firstTerm=sum((R-UV).*(-ones(N,1)*V(:,d)').*mask,2);
        secTerm=0.5*pDeriv'*denomCompU*U(:,d)*sigma_r^2;
        U_new(:,d)=U(:,d)-eta*(firstTerm+secTerm);
    end
    
    U=U_new;
    
    % Gradient descent on V;
     V_new = zeros(M, D);
     
     UV=U*V';
     pDeriv=eye(M);
     for d=1:D
         firstTerm=sum((R-UV).*(-U(:,d)*ones(1,M).*mask),1);
         secTerm=0.5*pDeriv'*denomCompV*V(:,d)*sigma_r^2;
         V_new(:,d)=V(:,d)-eta*(firstTerm'+secTerm);
     end
    
%     for n = 1 : N
%         for d= 1 : D
%             % Calculate the derivative w.r.t U(n,d)
%            
%             % First term
%             firstTerm = sum(mask(n,:).*(R(n,:) - U(n,:) * V').*(-V(:,d)'));
%             
%             % Second term
%             pDeriv = zeros(N,1);
%             pDeriv(n) = 1;
%             secTerm = 0.5 * pDeriv' * denomCompU * U(:,d) * sigma_r^2; % FIXME: sigma_r
%             U_new(n,d) = U(n,d) - eta * (firstTerm + secTerm);            
%         end
%     end

    
%     % Gradient descent on V
%     V_new = zeros(M, D);
%     for m = 1 : M
%         for d = 1 : D
%             % Calculate the derivate w.r.t V(m,d)
%             
%             % First term
%             firstTerm = sum(mask(:,m).*(R(:,m) - U * V(m,:)').*(-U(:,d)));
%             
%             % Second term
%             pDeriv = zeros(M,1);
%             pDeriv(m) = 1;
%             secTerm = 0.5 * pDeriv' * denomCompV * V(:,d) * sigma_r^2; % FIXME: sigma_r
%             V_new(m,d) = V(m,d) - eta * (firstTerm + secTerm);
%         end
%     end
    
    % Update U and V
    U = U_new;
    V = V_new;
    numIters = numIters + 1;
    
    % Calculate the error with current U and V
    tmpMat = mask.*((R - U*V').^2);
    err = 0.5 * sum(tmpMat(:)) / (sigma_r^2);
    for d = 1 : D
        err = err + 0.5 * (U(:,d)'*K_u_inv*U(:,d) + V(:,d)'*K_v_inv*V(:,d)) * sigma_r^2; % FIXME: sigma_r
    end
    E(numIters) = err;
    err

    if (err < 4e4)
        break;
    end
%     if (numIters > 1)
%         if (E(numIters) - E(numIters - 1) > 0)
%             disp('Error increases again...');
%         end
%         % Converges when the error stops changing for 4 consecutive times
%         if (((E(numIters -1) - E(numIters))/E(numIters - 1) < thresPercent && numIters > minIters) || numIters > maxIters)
%             if (numIters > maxIters)
%                 break
%             end
%             stopUpdateAtIter2 = numIters;
%             if (stopUpdateAtIter2 - stopUpdateAtIter1 == 1)
%                 stopTimes = stopTimes + 1; 
%                 if (stopTimes == 4)
%                     break
%                 end
%             else
%                 stopTimes = 0;
%             end
%             stopUpdateAtIter1 = stopUpdateAtIter2;
%         end
%     end
end